Midgap states in corrugated graphene: Ab-initio calculations and effective field theory 
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We investigate the electronic properties of corrugated graphene and show how rippling-induced 
pseudomagnetic fields alter graphene's low-energy electronic properties by combining first principles 
calculations with an effective field theory. The formation of flat bands near the Fermi level corre- 
sponding to pseudo-Landau levels is studied as a function of the rippling parameters. Quenched 
and relaxed ripples turn out to be fundamentally different is this respect: It is demonstrated, both 
numerically and analytically, that annealing of quenched ripples can destroy the fiat bands. 



Graphene, i.e. a monolayer of graphite, is the first 
truly two dimensional (2D) material available for ex- 
periments [H, [3] ■ Its low energy electronic structure re- 
sembling 2D Dirac massless fermion dynamics with the 
speed of light c being replaced by the Fermi velocity 
Vf ~ c/300 makes ultrarelativistic physics observable in 
this material 0, |3i Si- These extraordinary electronic 
properties are immediately related to the graphene's 2D 
honeycomb crystal structure, which was puzzling itself: 
Thermal fluctuations in two-dimensional solids, in prin- 
ciple, should lead to huge displacements of the carbon 
atoms from their perfect lattice arrangement and de- 
stroy any long range crystalline order. As an important 
step towards solving this puzzle, transmission electron 
microscopy ^] and atomistic simulations [7] found that 
free standing graphene sheets are not perfectly flat but 
exhibit ripples. In addition to this intrinsic crumpling, 
graphene's bonding to a substrate can also introduce rip- 
pling [1, Both, the intrinsic and extrinsic ripples, may 
be accompanied by different strengths of in plane strain. 

The effect of curvature and strain on the electronic 
structure of graphene can be described by effective gauge 
flelds A acting on the electrons 13, lH - in close ana logy 
to curvature effects in carbon nanotubes (CNTs) [H, [l3| . 
There, curvature decisively determines the low energy 
electronic properties, as it can open band gaps: The in- 
duced effective gauge fields shift the Fermi points of the 
2D graphene bands away from the ID nanotube bands. 
In graphene, on the contrary, the Dirac points are simi- 
larly shifted by uniform gauge flelds but no gap opening 
is expected unless a modulated electrostatic potential is 
present Rather, the Dirac point shiftin g le ads to 



unusually strong electron-phonon interactions }15| . 

However, the rippling induced gauge fields A are 
nonuniform and affect the electrons in g raphene like an 
effective mag netic field S = V x A [lO|. Tight-binding 
(TB) estimations on the effective magnetic field induced 
to graphene by rippling found the possibility of partially 
flat bands, which are the analog of Landau levels in real 



magnetic flelds Therefore, these flat bands are re- 

ferred to as pseudo Landau levels. In particular, zero- 
energy chiral states (n = LL) at the Fermi level should 
occur in inhomogeneous (real and effective) magnetic 
flelds. as was pointed out from topological considerations 
0i ■ The appearance of ripplc-induccd mid-gap states 
can lead to important consequences: The increased DOS 
at the Fermi level will enhance the tendency to spatial 
inhomogeneities |14| as well as lead to strong resonant 
electron scattering!^ 17, [l^. Knowing the conditions 
under which midgap states occur is therefore essential to 
resolve the debate on electron scattering in graphene. 

So far, the predictions on midgap states have been 
rather qualitative involving adjustable parameters, like 
hopping matrix elements, their change with strain and 
curvature as well as they neglected rehybridization ef- 
fects of the TT and a bands. It is not clear a priori how 
essential these effects can be. For example, taking into 
account next-nearest-neighbor hopping leads to a scalar 
electrostatic potential induced by the ripples [3] which 
can cause opening a gap around the Fermi level [l^ . 
Only based on complete flrst-principles calculations one 
can judge how important the corrections to the simplest 
nearest-neighbor TB model are. 

In this letter, we present full potential density func- 
tional theory (DFT) studies of quenched and annealed 
graphene ripples. We flnd flat bands very close to the 
Dirac point for quenched ripples, whereas in annealed 
ripples these midgap states turn out to be suppressed. 
Based on a nearest neighbor TB model, we extend the 
low energy effective fleld theory description of graphene 
to include this relaxation effect. The qualitative agree- 
ment of our effective fleld theory with DFT justifles this 
nearest neighbor TB based theory for describing corru- 
gated graphene. 

In our DFT calculations, quenched ripples are mod- 
elled as sinusoidal graphene ripples with height fleld 
h{x,y) — hosin{qx) (Fig. [U upper panel). While the 
effective gauge fleld A is caused by strain and curvature. 
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Figure 1: (Colour online) Upper panel: Schematic top and 
side view of the sinusoidal graphene ripples. The rippling pe- 
riod is denoted by I and their amplitude by h. Lower panel: 
Perspective view of a relaxed graphene ripple. The green 
arrows show the displacement of the atoms during the relax- 
ation. To enforce a constant rippling amplitude the vertical 
position of the atoms marked as big red dots has been fixed. 



the strain in the ripple can be significantly reduced by 
allowing all atoms and the supercell shape to relax with 
the constraint of fixed rippling amplitude Hq. As starting 
point for the relaxed cell shape, we shortened the super- 
cell in rippling direction such that the arc length of one 
sinusoidal graphene ripple period coincides with the equi- 
librium length of the same supercell for flat graphene. 
Then, standard relaxation of the atomic positions and 
the cell shape under the rippling constraint from above 
has been performed. This relaxation leads to our model 
of annealed graphene ripples, where "annealed" means 
that in these ripples possible external sources for energy 
barriers preventing the ripples from relaxation have been 
removed. Such sources might be im pur ities or bonding 
to a substrate with lattice mismatch. [20l [2l[ 

In all of our calculations we emp loy the generalized 
gradient approximation (GGA) [2^ to DFT for super- 
cells containing up to 160 carbon atoms. The resulting 
Kohn-Sham problems are solved with the Vienna Ab Ini- 
tio Simulation (VASP) [1^ package by expanding the 
electronic bands into projector augmented waves (PAW) 
24 , 2^ . Plane wave cut-offs of 500 eV for band-structure 
and 875 eV for the relaxations and total energy calcula- 
tions were used. For the total energy calculations, the 
Brillouin zone integrations were performed with 0.1 eV 
Gaussian smearing on k-meshes denser than 24 x 24 when 
folded back to the graphene first Brillouin zone, whereas 
for relaxations and input charge densities for band struc- 
ture calculations k-meshes diluted by a factor of 2 turned 
out to be sufficient. 

Firstly, the occurrence of the n = LL in quenched 
ripples as function of the ripple length / and ho/l has 
been studied, in this way. For ho/l = lA/4bo with 
bo = \/3ao/2 and oq = 2.465A being the graphene lat- 
tice constant, prominent changes in the high valance and 



Figure 2: (Colour online) Band structure of corrugated 
graphene sheets: The highest valence and lowest conduction 
bands along the ky direction (perpendicular to the rippling 
direction) are shown in main panel for different armchair rip- 
ples. Dash dotted and solid lines: Unrelaxed sinusoidal ripples 
with fixed ho/l — lA/4bo ratio for I = 16, SOfeo- Dashed lines: 
A ripple relaxed with the constraint ho = 4A and a rippling 
period of 16 graphene unit cells. The energy required to cre- 
ate the quenched ripples considered, here, is 0.5 eV/atom. 
For the relaxed ripple this energy is 0.03 eV/atom 



the low conduction bands occur (Fig. [5]) depending on 
the rippling period I. The shorter ripple {I — 166o) 
exhibits electron dispersion resembling massless parti- 
cles with the Dirac point shifted from ky = 27r/3ao to 
ky « 0.375(27r/ao), where ky is the crystal momentum 
perpendicular to the rippling direction, x. However, the 
bands of the I = 80 bo ripple are flat near the Fermi level 
and exhibit all characteristics of the n — LL of Dirac 
fermions: They are chiral, that is, fully sublattice polar- 
ized and localized in regions of maximum effective mag- 
netic field B. 

This manifests in the local density of states (LDOS) 
(see Fig. [3]) as follows. For a ripple of the form 
h(x,y) — hosiii{qx) with q — 2tt/1 and I sufficiently 
large, the effective gauge field is ^ ~ {hqcos{qx))^ and 
accordingly B ~ q'^ sm{2qx) . Thus, for I — 806o the 
absolute value of the effective field is maximum around 
X = 10, 30, 50, 706o. In this high field region (see Fig. [3]), 
the spectrum is gapped around the Fermi level (E = 0) 
in one sublattice (here B) but exhibits a mid-gap peak in 
the other sublattice. This will lead to sublattice stripes in 
low bias STM topographs with sublattice A bright and B 
dark, which is very reminiscent of midgap impurity states 
[i^t . In the low effective field region, the LDOS recovers 
the pseudogap shape typical for fiat graphene and the 
two sublattices appear fully equivalent again. 

The band structures of the quenched ripples shown 
in Fig. [2] allow to estimate the strength of the in- 
volved pseudomagnetic fields: Given the shift of the 
Dirac point of Sky — 0.042(27r/ao) away from the fiat 
graphene value, we obtain the average gauge field = 
h6ky/e — 2.8 • lO'^TflQ. Due to the sinusoidal shape 
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Figure 3: (Colour online) The local density of states (LDOS) 
inside the cells at x = 16o (low efF. field) and at a; = lOfeo 
(high field region). For the low field region, the LDOS is the 
same in both sublattices (only sublattice A plotted is here, 
dashed line), whereas in the high field region the LDOS in 
sublattice A (solid) and B (dash-dotted) differ significantly. 



of the ripple the amplitude sinusoidal pseudomagnetic 
field is Bo = Aqq w 250T for q = 27r/80&o. Due to 
B = Bo sin(2(7a;) a pseudo landau level wave function 
should be localized on a length less than I /A correspond- 
ing to an area of /^/16 = 18 nm^. (See also Ref. 
250T correspond to approx. 2 (pseudo)magnetic flux 
quanta per ISnin^. 

The ripple under consideration has maximum local 
strain of 24%, which is more than the avera ge s train of 
approx. 1% measured in epitaxial graphene 
found in other nanosized epitaxial materials 
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or 4% 
Scal- 
ing down the pseudo magnetic field to these strain values 
yields lOT and 40T, respectively. 

So far, we considered quenched ripples, but the pseu- 
domagnetic field is sensitive not only to flexural defor- 
mations but also to in-plane distortions. In general, rip- 
ples will be accompanied by in-plane distortions: For the 
quenched ripple of length I — 165o with h/l ~ lA/4bo, 
we relaxed the atomic positions and the supercell shape 
of the ripple with the constraint of fixed rippling height. 
This "annealed" ripple is still sinusoidal (Fig. [1] lower 
panel) but with all nearest neighbor bond lengths be- 
ing equal after relaxation. During this relaxation process 
the effective gauge field decreases as can be seen from 
the band structure of the relaxed ripple in Fig. [2l The 
Dirac point for the relaxed ripple is at ky w 0.332(27r/ao), 
which is by a factor of 40 closer the flat graphene value 
of ky = l/3(27r/ao) than the Dirac point of the quenched 
ripples. This corresponds to a decrease of the average 
effective gauge field the same factor. This suppression of 
effective gauge and pseudomagentic fields in the annealed 
ripple can be understood in terms of the following model: 
Consider the nearest neighbor TB Hamiltonian 



tnn{k) 



(1) 



the Fermi operators of electrons in sublattice A (B) with 
crystal momentum k. Close to ±K = (0, =F47r/3-\/3), in 
the corner of the Brillouin zone, the dispersion for the un- 
deformed lattice vanishes linearly and a continuum the- 
ory describing low energy electronic states can be defined. 
As in Ref. [28|, we introduce a pair of two-dimensional 

spinors = ^ ^ ' ^2 = ^ ^describing 

electronic wave packets centered at K and — K point, 
respectively, and expand the hopping integrals. Near the 
K point the undeformed hopping integral is 



tnn (k) 



3toao 



{qy - iqx) 2vo5, 



(2) 



where q — k~K, B — ■^{dx + idy) and vo = 3toS.o/2 is the 
Fermi velocity, involving explicitly the nearest-neighbor 
spacing ao = ao/^/3. Slow lattice deformations with the 
deformation tensor field 



Uab = -^[daUb + dbUa] + dahdbh 



(3) 



with a,b Cz {x, y} lead to changes in the hopping integral 



where ei = (1,0), ea = (1/2, \/3/2), eg = (1/2,-^3/2) 
are vectors connecting nearest neighbor sites and z = 
x+iy. Thus the nearest neighbor hopping matrix element 
near the point K is now 



2vo {B - 



(5) 



where 7 — dlnto/dR plays a role of the charge, and Uzz 
can be interpreted as a vector potential: 



A"^ — Uzz: A"^ — A'l = A't iA?. — Uzz 



(6) 



where t„„ are the hopping matrix elements and Ak (Bk) 



We have considered smooth lattice deformations and 
established that they generate an Abelian vector poten- 
tial having opposite signs for different valleys. It is worth 
noticing that this vector potential constitutes a part of 
the non- Abelian field whose other noncommuting compo- 
nents A^'^ are generated by abrupt changes of the near- 
est neighbor hopping integrals. The complete low energy 
Hamiltonian is |28| 

H = {l<S)IAo + voa^ ® [-i% + 7t"A;^]} * (7) 

where the Pauli matrices (t^,/! — x,y act in the spinor 
space and matrices r", a = 1, 2, 3 act on the valley index 
of the spinor \E'+ = {'^^,^^). 

The zero energy wave functions for the vector potential 
deformations can be found analytically for Ao = 0. For 
simplicity we will restrict ourselves to the case of slow 
deformations when A^'^ = 0. Then the zero energy wave 
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functions (if they exist) are expressed in terms of pseudo- 
magnetic field 



For B > we have: 



A 



K 



K 



z exp 



dd 



B 



Bk = S_ 



K 



(9) 



and for B < Q one has to interchange Ak,A-k with 
Bk,B^k and z with z. The power n ranges from to 
the integer part of the pseudomagnetic flux. As many 
zero modes exist. 

The expression for B and for the deformation tensor 
changes drasticahy depending on whether the elastic en- 
ergy is at its minimum or not. The expression for the 
elastic energy density of a smooth surface compatible 
with the C3 symmetry is given by [29| 



E 



(A + m) 



{du + du) + dhdh 



K 



+fi [du + {dhf] [du + {dh)^] + —{V^h) 



^iAA + (A + ^) 



(10) 



where TZ — [d^hd^h—{ddh)^] is the Gaussian curvature of 
the surface. If the membrane is at equilibrium, the elastic 
energy (|10p is at its minimum and the vector potential is 
given by 



A^ 



(A + /i) 
(A + 2fi) dd 



n. 



(11) 



TZ is the Jacobian of the coordinate transformation ^ 
dh, l^dh: 



d^hd^h - [ddhf 



d{z,z) 



(12) 



Thus, TZ and, as a consequence, A'^^A^ as well as ^0 
vanish on any configuration of h which depends on just 
one Cartesian coordinate. This includes all plane wave 
configurations, in particular those studied by DFT in this 
letter. 

Substituting Eqn. (fTT|) into ([5]) we get 



B 



i(A + iJL)d^- d^ 
X + 2fi {ddy 



[d^hd^h - {ddhf 



(13) 



which is C3 symmetric, as it must be. This C3 sym- 
metry of B for a relaxed membrane leads to important 
qualitative differences with a real magnetic field. A local 
vortex of magnetic field with flux TV carries N normaliz- 
able zero modes, but to create such a flux by deforming 



a membrane is not possible. An analogue of a magnetic 
vortex is a point-like defect h = /i(|r — a|) carrying lo- 
cal Gaussian curvature TZ — TZ{\r — a|). Then the argu- 



(8) ment in wave function ^ is -^B — 7 sin(3(/))6(| 



where b{r) = J (2^)fc^ J-i{kr)TZk and (/> is the angle with 
respect to the crystalline axis. Since TZk is constant at 
small fc, at distances larger than the size of the defect 
6(r) ~ r. Such wave function is not normalizable. On 
the other hand the wave function with two defects with 
curvature of the opposite sign will be normalizable for 
the entire volume (as a plane wave). Such state is non- 
degenerate, that is we have n = in So, if we allow 
the graphene membrane to relax there are only nonde- 
generate zero modes, which only exist for membrane con- 
figurations with nonzero Gaussian curvature. Therefore, 
the effect of relaxations on the electronic properties of 
the graphene ripples is qualitatively the same for ID and 
2D ripples. The degenerate zero modes are suppressed 
in relaxed ripples and no significant midgap peak in the 
total DOS of 2D ripples is expected to occur. 

In conclusion, our results demonstrate an essential dif- 
ference between quenched and annealed ripple structure. 
If the system is allowed to relax to its minimum of elas- 
tic energy for a given h{x,y) profile it decreases dras- 
tically the amplitude of pseudo-magnetic fields and can 
lead to disappearance of the mid-gap states. It may be 
an important statement in light of the hypothesis [l3| 
that quenched ripples are the main source of scattering in 
graphene and the very recent observation that annealing 
of a freely hanged grap hene membrane can increase dras- 
tically its mobility [30| : Upon annealing, impurities caus- 
ing energy barriers that prevent the suspended graphene 
from relaxation might be removed. This issue requires 
further investigations. Our ab-initio calculations reveal 
a perfect particle-hole symmetry at low energies and jus- 
tify the nearest neighbor hopping based field theory for 
describing the electronic properties of graphene ripples. 
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